int power(int base, int exponent)
{
  int r = 1;
  int t;
  for(t = exponent ; t>0 ; t--)
    r = r*base;
  return r;
}

void CubeStats (int *first, int edgeLength, int dimension, int bArraySize,
                int *saveRange, int *saveAverage) {

    /* Array of coefficients, initially all 0s. */
    int ArrayOfCs[dimension];
    int index;
    for (index = 0; index < dimension; ++index) {
        ArrayOfCs[index] = 0;  }

    int i, j, k, q;
    int subCubeRange;

    int subCubeAverage;
    int * polynomial;
    for (i = 0; i < power(edgeLength, dimension-1); ++i) {
        for (j = 0; j < edgeLength; ++j) {
            ArrayOfCs[0] = j;       /* Increment coefficient 1. */

            /* Polynomial := initial position, add one term of the polynomial at a time ( c_{n} b^{ d_{n}-1 } ). */
            polynomial = first;
            for (k = 0; k < dimension; ++k) {
                polynomial += ArrayOfCs[k] * power(bArraySize, k);  }

            if (*polynomial < min) {
                min = *polynomial;  }
            else if (*polynomial > max) {
                max = *polynomial;  }
            total += *polynomial;  }

        /* Update coefficients. */
        q = 1;
        while (ArrayOfCs[q] == edgeLength-1) {
            ArrayOfCs[q] = 0;
            q++;  }
        ArrayOfCs[q]++;  }

    subCubeAverage = (int) (total / power(edgeLength, dimension));
    *saveAverage = subCubeAverage;

    subCubeRange = max - min;
    *saveRange = subCubeRange;
}
